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Abstract 

This paper is concerned with the problem of finding a quadratic common Lyapunov function for 
a family of stable linear systems. We present gradient iteration algorithms which give deterministic 
convergence for finite system families and probabilistic convergence for infinite families. 
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1 Introduction 

A switched system is a dynamical system described by a family of continuous-time subsystems and a rule 
that governs the switching between them. Such systems arise as models of processes regulated by switching 
control mechanisms and/or affected by abrupt changes in the dynamics. It is well known and easy to 
demonstrate that switching between stable subsystems may lead to instability. This makes stability analysis 
of switched systems an important and challenging problem, which has received considerable attention in 



the recent literature. We refer the reader to the survey articles [12] and || for an overview of available 
results. 

Stability of a switched system is guaranteed if the individual subsystems share a common Lyapunov 
function. If the system enjoys this property, then stability is preserved for arbitrary switching sequences. 
When the subsystems being switched are obtained as feedback interconnections of a given process with 
different stabilizing controllers, this means that one does not need to worry about stability and can con- 
centrate on other issues such as performance. 

A particular case of interest is when the subsystems are linear time-invariant and a quadratic common 
Lyapunov function is sought. Although a number of conditions for the existence of such a Lyapunov function 
have been obtained, general results are lacking. On the other hand, the problem of finding a quadratic 
common Lyapunov function amounts to solving a system of linear matrix inequalities (LMIs) , and efficient 
methods for solving such inequalities are available [||. However, this approach becomes infeasible when 
the number of subsystems being switched is large, and is not useful in the case of an infinite family of 
subsystems. 
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In view of the above remarks, there is a need for developing computationally tractable algorithms 
which can be used to find quadratic common Lyapunov functions for large families of linear systems. It 
turns out that this can be achieved by handling matrix inequality constraints sequentially rather than 
simultaneously (For comparison, we mention the iterative algorithm proposed in [^] for a special case of 
the above problem, namely, finding a common Lyapunov function for a finite family of commuting linear 
systems.) In the case of an infinite family of systems, one can relax the objective and apply the same ideas 
combined with randomization to obtain probabilistic convergence results. 

The use of gradient algorithms for solving matrix inequalities was proposed in P] and in the 
context of probabilistic control design for uncertain linear systems; see also g g H, II for related 
recent developments. The main ideas go back to the early work on solving algebraic inequalities reported 
in |l4|, |l7|, 20|. The goal of this paper is to show how the problem posed above can be addressed 



using similar techniques. We describe a general framework for developing gradient iteration algorithms 
for finding quadratic common Lyapunov functions and discuss and compare several specific design choices, 
some of which were used in previous work cited above while others were not. Convergence to a quadratic 
common Lyapunov function (if one exists) in a finite number of steps is guaranteed for a finite family of 
linear systems and takes place with probability one for an infinite family. While our approach parallels 
that of [p.8y , there are several important technical differences on which we elaborate below. We first treat 
finite families and then discuss an extension to infinite families. 



2 Problem formulation and notation 

Suppose first that we are given a family of real Hurwitz n x n matrices A\, . . . , An, where n and N are 
positive integers. We write P > (or P > 0) to indicate that a matrix P is symmetric nonnegative definite 
(respectively, positive definite), and P < (or P < 0) to indicate that P is symmetric nonpositive definite 
(respectively, negative definite). We assume throughout that there exists a matrix P > which satisfies 

PAt + AjP <0, i = l,...,N. 

This means that the quadratic function V(x) := x T Px is a common Lyapunov function for the family of 
asymptotically stable linear systems 

x = AiX, i = l,...,N. (1) 

Fix an arbitrary matrix Q > 0. Multiplying P by a sufficiently large positive number, we see that the 
system of inequalities 

PAi + AfP + Q<0, i = l,...,N (2) 

has a solution P > 0. Moreover, if a symmetric matrix P satisfies (at least one of) the inequalities (|2|), 
then it is well known that we automatically have P > 0; see, e.g., p. 132]. Thus the problem of finding 
a quadratic common Lyapunov function for the family ([]]) is equivalent to that of finding a symmetric 
matrix P which satisfies (|2|). In what follows, we will be concerned with the latter problem. We denote 
the set of symmetric solutions of the inequalities @ by £. 

The space of real symmetric nxn matrices is equipped with the Frobenius norm | |J2| | := (Y17j=i R-ij) 1 ^ 2 
and the inner product (R, S) := tr RS. It is important to note that since the set C is nonempty by 
assumption, it must have a nonempty interior. Indeed, if P G C, then a standard perturbation argument 
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such as the one given in || Example 5.1] can be used to show that C contains a neighborhood of 7P for 
every 7 > 1. In fact, we have 7P + aP € £ for every symmetric matrix aP satisfying the bound 

A max (AP)< (T-l)WQ) 



max^o- max (Ai 



where A max and A m ; n denote the largest and the smallest eigenvalue, respectively, and <7 max denotes the 
largest singular value. It is well known that A max (AP) < ||aP||; see, e.g., @, pp. 296-297]. Since 7 can be 
arbitrarily large, we see that for every r > there exists a ball of radius r which is contained in C. This is 



in contrast with the case of Riccati inequalities treated in [18|, where a sufficiently small r with the above 
property is assumed to exist but is not explicitly known. 

We will need the notion of projection of a symmetric matrix R onto the convex cone of nonnegative 
definite matrices. This projection is defined as 

R + := argmin \ \R — S\ I. 

S>0 



It can be computed explicitly as follows (see |l7|| ): If R = UAU T , where U is orthogonal and A is diagonal 
with entries Ai, . . . , A n , then R + = UA + U T , where A + is diagonal with entries max{0, Ai}, . . . , max{0, A n }. 
We also denote R — R + by R~ . 



3 Gradient algorithms 

Let / be a convex functional on the space of symmetric matrices, which assigns to a matrix R a real number 
f(R), with the property that f(R) < if and only if R < 0. We suppose that this functional is differentiable 
(this condition can be relaxed somewhat, as discussed later) and that its gradient is symmetric. We denote 
this gradient by dnf, and similar notation will be used for other gradients appearing below. Specific 
examples of functionals with the above properties will be given in Section [|; here we keep the discussion 
general. 

Given a symmetric matrix P and another matrix A, we let 

v(P,A) := f(PA + A T P + Q) 

where / is the functional introduced above and Q > is the matrix from (|2|). Since / is convex, v is 
convex in P. Well-known results imply that for each integer i between 1 and N, solutions of the gradient 
system P = -d P v(P, A { ) converge to the set {P : v(P, Ai) < 0} = {P : PA { + AjP + Q < 0}. The same is 
true for the associated discrete iterations Pk+i = Pk — fJ-kdpv(Pk, Ai), k = 0, 1, . . . if the step-sizes are 
chosen appropriately. Moreover, by switching between the above iterations for different values of i, we can 
make Pk converge to the intersection of the corresponding sets, which is precisely the set C of solutions 
of (|2|). This happens because the distance from P^ to £ with respect to the Frobenius norm is a decreasing 
function of k. 

We now make the above discussion precise by describing how the gradient iterations are to be carried 
out to ensure convergence to C in a finite number of steps. We need to pick a "scheduling function" h from 
nonnegative integers to the set {1, . . . , N} which has the following revisitation property: For every integer 
i between 1 and N and for every integer I > there exists an integer k > I such that h{k) = i. An obvious 
choice is h(k) := (k mod N) + 1 = k - N[k/N\ + 1. 
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For each k = 0, 1, . . . , let us define the step-size by the formula 



_ av{Pk, A h{k) ) + r\\d P v(P k ,A h{k) )\\ 

\\d P v(P k ,A h{k) W [ ' 

where < a < 1 and r > are arbitrary. Consider the iterations 

P k+1 = S Pk ~ ^dpv(P k ,A h(k) ) if v(P k ,A h{k) ) > 0, 
1 P k otherwise. 

We take the initial condition Pq to be symmetric. (For example, a solution of one of the inequalities @ 
provides a convenient choice for Pq.) Then P k is symmetric for each k, since dpv is symmetric in light of 
the following lemma and the fact that d R f is symmetric. 

Lemma 1 The gradient of v is given by 

d P v(P, A) = Ad R f(PA + A T P + Q)+ 8 R f(PA + A T P + Q)A T . (5) 

Proof. Denoting by aP a small perturbation in P and by ~ equality up to first-order terms in aP, we 
write 

v(P + aP, A) = f(PA + A T P + Q + aPA + A t aP) 

» f{PA + A T P + Q) + {d R f(PA + A T P + Q),aPA + A t aP) 
= v(P, A) + tr [[Ad R f{PA + A T P + Q) + d R f(PA + A T P + Q)A t )aP 

from which (|5|) follows. □ 

On the other hand, P k +i is not guaranteed to be positive definite or at least nonnegative definite, even 
if P k > 0. To make sure that a nonnegative definite matrix is generated at every step, we could instead 
consider 

P k+1 = [^ Pk ~ Vkdpv(P k ,A h ( k) )} + if v(P k , A h(k) ) > 0, 
I P k otherwise 

using the projection operation defined in Section pi. Since all matrices in the set C are positive definite, 
this modification also has the potential of improving convergence; see the proof of Theorem [l] below. 

Remark 1 The algorithm (||) exactly parallels the one proposed in for finding nonnegative definite 
solutions of systems of Riccati inequalities. The Lyapunov inequalities (Q), on the other hand, have the 
property that if a symmetric matrix P k satisfies at least one of them, then necessarily P k > 0. This means 
that the projection is not really needed, and the convergence result presented in the next section implies 
that the algorithm (|j) generates only nonnegative definite matrices after sufficiently many steps. 



4 Deterministic convergence for finite families 

We now demonstrate that the above gradient algorithms indeed provide convergence to the desired set L 
in a finite number of steps. When for a given k we have v(P k , A^^) > 0, we say that a correction step is 
executed. 
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Theorem 1 When the algorithm (0) or the algorithm (|6|) is applied with the step-size given by (^), there 
exists an integer k* such that P k * € C 



Proof. The proof is carried out along the lines of [18]. Consider a k for which v(P k , A-hVk)) > and so a 
correction step is executed. As shown in Section [|, the set C contains a ball of radius r, centered at some 
matrix P*. We will prove that 

||P fc+ i-Pl 2 < ||P fc -P*|| 2 -r 2 . (7) 

Since the revisitation property of h guarantees that a correction step occurs at least once in every N steps 
until P k G C, we can conclude that no more than iV[J[Po — P*|| 2 /V 2 J steps are needed, and the proof will 
be complete. 

Consider the algorithm (|6j). We have 

||P fe+1 - P*|| 2 = \\[P k - ti k d P v(P k ,A m )]+ - P*\\ 2 < \\P k - fi k d P v(P k ,A m ) - P*\\ 2 

where the last inequality follows from the definition of the projection and the supporting hyperplane 
theorem. Thus we see that it is enough to show ([/]) for the algorithm (|j). 

To this end, define 

P := P * + iTa (J a Tu d P v ( P k> A h(k)) e C - 

Then we write 

||P fc+1 - P*|| 2 = ||P fc - P* - v k d P v(P k ,A m )\\ 2 = \\P k - P*\\ 2 + n 2 \\d P v(P k ,A m )\\ 2 
- 2» k (d P v(P k ,A m ),P k - P) - 2fi k {dpv(P k ,A m ),P- P*). 

We now consider the last two terms. Due to convexity of v in P, we have 

(d P v(P k , A h(k) ), P k - P) > v(P k , A h{k) ) > av{P k ,A h(k) ) 

while the definition of P gives 

(d P v(P k , A h(k) ),P-P*)=r\\dpv(P k ,A h(k) )\\. 

Therefore, 

||P fe+ i-P*|| 2 < ||P fc -Pl| 2 + ^ 2 ||^(P fc ,A Mfc) )|| 2 -2^(a U (P fe ,A h(fc) )+r||9pKPfe,^ 
Substituting the value of fi k defined in (||), we obtain 

HP P*||2/||p p*||2 (^(^,^(fc))+^ll^(^^Mfc))H) 2 ^ up p*||2 2 

\\Pk+i ~ P \\ < \\Pk ~ P \\ — — z run S\\P k -P\\ -r 

\\dpv(P k ,A h{k) )\\ 2 

and so the inequality (^) holds as claimed. □ 

Since the matrix P* used in the above proof is not known, the number k* may be difficult to estimate 
in practice. One can of course let the algorithm run for some number of steps k and then check whether 
or not the matrix P k satisfies the inequalities (g); performing such a check is an easy task compared with 
that of solving the inequalities directly. We also have a lot of freedom in the choice of the step-size; more 
information on how to choose the step-size efficiently in gradient algorithms of this kind can be found, e.g., 
in 0,0,^. 
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5 Possible choices for / 



We suggest two admissible choices of the functional /. One is given by 

f(R):=\\R + \\ 2 . (8) 

This functional has all the properties required in Section [3|, as shown in jlT]]. For completeness, we sketch 
the argument below. 

Lemma 2 The functional (|8|) is convex and differentiate, with gradient given by 

d R f(R) = 2R+. (9) 

Proof. The calculations that follow rely on some standard properties of the projection and the cone of 
symmetric nonnegative definite matrices; see |17j for details. We write 

f(R + AR) = \\R + AR-(R + aR)~ II 2 = inf \\R + aR- S\\ 2 

S<0 

< \\R + aR- R-\\ 2 = \\R + +AR\\ 2 = \\R + \\ 2 + {2R + ,aR} + \\aR\\ 2 
« f(R) + (2R + ,AR). 

On the other hand, 

f(R + aR) = \\R + aR-(R + aR)- || 2 = \\R + + R~ + aR - (R + aR)~ \\ 2 

= \\R + \\ 2 + (2R + ,aR) + 2(R + ,R-} + 2{R + ,-(R + aR)~) + \\R~ + aR - (R + AR)~ \\ 2 
> f(R) + {2R + ,AR). 

This proves (||), and the inequality f(R + aR) > f(R) + (djif(R),AR) implies that / is convex. □ 
Combining the formulas @ and (P), we see that with this choice of / the gradient of v is given by 

d P v{P, A) = 2A{PA + A T P + Q)+ + 2{PA + A T P + Q) + A T . 

This is the quantity that needs to be calculated at every step when implementing the algorithms of Section ||. 
As explained in Section ^, this is a routine calculation based on eigenvalue decomposition. An analogous 
construction was used in Ifii 



Another option is to let f(R) be the largest eigenvalue of R: 

f(R) := A max (ii). (10) 
The following result is standard; see, e.g, @, p. 372]. 



Lemma 3 The functional ( |i0[) is convex and, when X maK (R) is a simple eigenvalue, differentiate with 
gradient given by 

d R f(R) = xx T 

where x is a unit eigenvector of R with eigenvalue A max (-R). 

With this choice of /, the gradient of v is given by 

d P v(P, A) = Axx T + xx T A T 

and the projection is no longer required for implementation of the algorithm (|j). A disadvantage of this 
/ is that it is not differentiable when X max (R) is not a simple eigenvalue. However, since a generic matrix 
has distinct eigenvalues, this problem can always be fixed by slightly perturbing the matrix P^ if necessary. 
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6 Extension: probabilistic convergence for infinite families 



Let us now consider the more general situation where we are given a compact set of real Hurwitz matrices 
A := {A p : p G V}, parameterized by some index set V which is in general infinite. It is still true that the 
problem of finding a matrix P > which satisfies the inequalities 

PA P + A T V P < Vp G V 

is equivalent to the problem of finding a symmetric matrix P which satisfies the inequalities 

PAp + A T V P + Q < Vp G V 

where Q > is arbitrary. As before, we denote the set of such matrices P by C and assume that it is 
nonempty; this again implies that C contains balls of arbitrarily large radii. 

When the set of matrices is infinite, choosing a matrix at each step using a function h with the 
revisitation property as in Section || is no longer possible. Instead, we use randomization. Namely, for 
each k we randomly pick a matrix in A according to some probability distribution on A with the following 
property: Every subset of A which is open relative to A has a nonzero probability measure. We let h(k) 
be the index of the matrix chosen in this way (this index may not be unique if the map p i— > A p is not 
injective). Note that A may be composed of several disjoint subsets such as intervals or isolated points, 
and we must assign a positive probability measure to each of them. If A consists of a finite number of 
isolated points, then randomization provides an alternative to the approach described in Section ||. As is 
well known, for infinite families described by convex polyhedra the problem reduces to the corresponding 



problem for finite families generated by the vertices; see, e.g., [13| and the references therein. 

The randomized versions of the algorithms from Section || can now be defined by the same formulas (Q) 
and © as before. Remark |l| also applies here. If ^ C, then v(Pk,A) > for some A G A, and by 
continuity the same is true for all A in a sufficiently small neighborhood of A in A. Since the probability 
measure of this neighborhood is positive, a correction step will be executed with probability one after 
a finite number of steps. Proceeding exactly as in the proof of Theorem [l], we arrive at the following 
probabilistic counterpart of that theorem. 

Theorem 2 When the algorithm (Q) or the algorithm (JgJ) is applied with the step-size given by (^j, the 
probability of obtaining Pk* G C for some integer k* is equal to one. 



The above result parallels Theorem 1 from |18|| , although the context here is different. In [18] it is also 
shown how one can compute a lower bound on the probability of finding a solution in a given number of 
steps, provided that additional a priori information is available. 



7 Concluding remarks 

We presented gradient iteration algorithms for finding a quadratic common Lyapunov function for a family 
of asymptotically stable linear systems. We derived a deterministic convergence result for finite families, 
and then used randomization to obtain a probabilistic counterpart for infinite families. The algorithms 
were described on a general level with some freedom left in the design choices, several possibilities for which 
were also discussed. 

We implemented the above algorithms in MATLAB for interval families of 4 x 4 and 5x5 upper- 
triangular Hurwitz matrices, for which quadratic common Lyapunov functions are known to exist (see, 
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e.g., JTT|). We applied the method of Section pi noting that it suffices to work with the vertices. This 
was a straightforward task, and in all cases we observed convergence in a reasonable number of iterations: 
approximately 5,000 for the 4x4 case and 75,000 for the 5x5 case. The number of actual correction steps 
was an order of magnitude lower than the total number of iterations. Trade-offs among different choices 
of the specific gradient algorithm and the step-size remain to be understood. It is important to note that 
no analytical results on the existence of a common Lyapunov function are available for generic matrices of 
dimension higher than 2x2, and that interval families of 4 x 4 and 5x5 triangular matrices give rise to 
systems of 2 10 and 2 15 LMIs, respectively, solving which simultaneously is computationally intractable. 

The results reported here are valid under the assumption that a quadratic common Lyapunov function 
exists. They do not offer insight into the question of how to determine whether or not such a Lyapunov 
function can be found; this issue needs further investigation. By combining the above algorithms with 
a suitable recursive partitioning procedure, we can compute common Lyapunov functions for subsets of 
the given family of systems when the overall common Lyapunov function may not exist or is not found 
in a prescribed number of steps. Stability of the switched system is then guaranteed under appropriate 



constraints on the switching rate between different subsets; cf. [12]. Problems for future work also include 
choosing an optimal schedule for the iterations with respect to some performance criterion (for example, 
the speed of convergence) and incorporating additional constraints on the Lyapunov function to be found 
(such as the average or worst-case decay rate). 
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